﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using SlimDX;

namespace Blake.Fluid
{
    public class PointSet
    {
        #region Field

        protected Random random = new Random();

        protected List<FluidParticle> particles = new List<FluidParticle>();
        public List<FluidParticle> Particles
        {
            get
            {
                return particles;
            }
        }

        protected List<RectangleParticle> bodies = new List<RectangleParticle>();
        public List<RectangleParticle> Bodies
        {
            get
            {
                return bodies;
            }
        }

        protected int frame;

        protected float time;

        protected List<int> grid;
        protected List<int> gridCount;

        protected int currentParticle = -1;

        protected float gridCellSize;

        protected const int maxNeighbors = 80;

        protected ushort[] neighborCount = new ushort[ushort.MaxValue];
        protected ushort[][] neighbor = new ushort[ushort.MaxValue][];
        protected float[][] neighborDist = new float[ushort.MaxValue][];

        #endregion

        #region Properties

        public int NumPoints
        {
            get
            {
                return particles.Count;
            }
        }

        public float DeltaTime { get; protected set; }

        public Vector3 GridResolution { get; protected set; }

        public Vector3 GridSize { get; protected set; }

        public Vector3 GridMinimum { get; protected set; }

        public Vector3 GridMaximum { get; protected set; }

        public Vector3 GridDelta { get; protected set; }

        public FluidSystemParams Parameters { get; protected set; }

        #endregion

        #region Events

        #endregion

        #region Constructors

        public PointSet()
        {
            grid = new List<int>();
            gridCount = new List<int>();

            GridResolution = new Vector3(0, 0, 0);

            for (int i = 0; i < ushort.MaxValue; i++)
            {
                neighborCount[i] = 0;
                neighbor[i] = new ushort[maxNeighbors];
                neighborDist[i] = new float[maxNeighbors];
            }

            Reset();
        }

        #endregion

        #region Public Methods

        #endregion

        #region Virtual Methods

        public virtual void Initialize(int numParticles)
        {

        }

        public virtual void Draw(Matrix viewMatrix, float radius)
        {
            //Loop through points
            //Use point position and color to draw a sphere
        }

        public virtual void Reset()
        {
            time = 0;
            DeltaTime = 0.1f;
            if (this.Parameters == null)
                this.Parameters = new FluidSystemParams();
            Parameters.EmitterRate = new Vector3(0, 0, 0);
            Parameters.EmitterPosition = new Vector3(0, 0, 10);
            Parameters.EmitterAngle = new Vector3(90, 45, 50);
            Parameters.EmitterDAngle = new Vector3(0, 0, 0);
            Parameters.EmitterSpread = new Vector3(4, 4, 1);
        }

        public virtual int AddPoint()
        {
            int index = particles.Count;
            particles.Add(new FluidParticle());
            return index;
        }

        public virtual int AddPointReuse()
        {
            //TODO: If at max elements, return random existing point instead?
            return AddPoint();
        }

        public FluidParticle GetPoint(int index)
        {
            return particles[index];
        }

        public virtual float GetValue(float x, float y, float z)
        {
            float dx, dy, dz, dsq;
            float sum = 0;
            int pointIndex;
            FluidParticle currentPoint;
            float R2 = 1.8f * 1.8f;

            var gridCells = Grid_FindCells(new Vector3(x, y, z), gridCellSize / 2.0f);

            foreach (int cellValue in gridCells)
            {
                if (cellValue != -1)
                {
                    pointIndex = grid[cellValue];
                    while (pointIndex != -1)
                    {
                        currentPoint = particles[pointIndex];
                        dx = x - currentPoint.Position.X;
                        dy = y - currentPoint.Position.Y;
                        dz = z - currentPoint.Position.Z;
                        dsq = dx * dx + dy * dy + dz * dz;
                        if (dsq < R2)
                            sum += R2 / dsq;
                        pointIndex = currentPoint.Next;
                    }
                }
            }
            return sum;
        }

        public virtual Vector3 GetGradient(float x, float y, float z)
        {
            Vector3 norm = new Vector3(0, 0, 0);
            float dx, dy, dz, dsq;
            int pointIndex;
            FluidParticle currentPoint;
            float R2 = (gridCellSize / 2.0f) * (gridCellSize / 2.0f);

            var gridCells = Grid_FindCells(new Vector3(x, y, z), gridCellSize / 2.0f);

            foreach (int cellValue in gridCells)
            {
                if (cellValue != -1)
                {
                    pointIndex = grid[cellValue];
                    while (pointIndex != -1)
                    {
                        currentPoint = particles[pointIndex];
                        dx = x - currentPoint.Position.X;
                        dy = y - currentPoint.Position.Y;
                        dz = z - currentPoint.Position.Z;
                        dsq = dx * dx + dy * dy + dz * dz;
                        if (dsq > 0 && dsq < R2)
                        {
                            dsq = 2.0f * R2 / (dsq * dsq);
                            norm.X += dx * dsq;
                            norm.Y += dy * dsq;
                            norm.Z += dz * dsq;
                        }
                        pointIndex = currentPoint.Next;
                    }
                }
            }
            norm.Normalize();
            return norm;
        }

        public virtual Color3 GetColor(float x, float y, float z)
        {
            Vector3 color = new Vector3(0, 0, 0);
            float dx, dy, dz, dsq;
            int pointIndex;
            FluidParticle currentPoint;
            float R2 = (gridCellSize / 2.0f) * (gridCellSize / 2.0f);

            var gridCells = Grid_FindCells(new Vector3(x, y, z), gridCellSize / 2.0f);

            foreach (int cellValue in gridCells)
            {
                if (cellValue != -1)
                {
                    pointIndex = grid[cellValue];
                    while (pointIndex != -1)
                    {
                        currentPoint = particles[pointIndex];
                        dx = x - currentPoint.Position.X;
                        dy = y - currentPoint.Position.Y;
                        dz = z - currentPoint.Position.Z;
                        dsq = dx * dx + dy * dy + dz * dz;
                        if (dsq < R2)
                        {
                            dsq = 2.0f * R2 / (dsq * dsq);
                            color.X += currentPoint.Color.Red * dsq;
                            color.Y += currentPoint.Color.Green * dsq;
                            color.Z += currentPoint.Color.Blue * dsq;
                        }
                        pointIndex = currentPoint.Next;
                    }
                }
            }
            color.Normalize();
            return new Color3(color.X, color.Y, color.Z);
        }

        public virtual void Run()
        {
            if (Parameters.EmitterRate.X > 0 && ++frame >= (int)Parameters.EmitterRate.X)
            {
                frame = 0;
                Emit(1.0f);
            }
            Advance();
        }

        protected virtual void Advance()
        {
            Vector3 velocityNext = new Vector3();
            Vector3 acceleration;
            Vector3 norm = new Vector3();

            velocityNext = Parameters.EmitterDAngle;
            velocityNext *= DeltaTime;
            Parameters.EmitterAngle += velocityNext;

            foreach (var particle in particles)
            {
                acceleration = new Vector3(0, 0, 0);

                if (Parameters.PlaneGravityEnabled)
                {
                    acceleration += Parameters.PlaneGravityDirection;
                }
                if (Parameters.PointGravityEnabled)
                {
                    norm = particle.Position - Parameters.PointGravityPosition;
                    norm.Normalize();
                    norm *= Parameters.PointGravityStrength;
                    acceleration -= norm;
                }

                //Leapfrog integration
                velocityNext = acceleration;
                velocityNext *= DeltaTime;

                // v(t+1/2) = v(t-1/2) + a(t) dt
                velocityNext += particle.Velocity;

                particle.VelocityEvaluative = particle.Velocity;
                particle.VelocityEvaluative += velocityNext;

                // v(t+1) = [v(t-1/2) + v(t+1/2)] * 0.5		used to compute forces later
                particle.VelocityEvaluative *= 0.5f;

                particle.Velocity = velocityNext;
                velocityNext *= DeltaTime;

                // p(t+1) = p(t) + v(t+1/2) dt
                particle.Position += velocityNext;
            }

            time += DeltaTime;
        }

        public virtual void Emit(float spacing)
        {
            FluidParticle particle;
            Vector3 direction = new Vector3();
            Vector3 position = new Vector3();
            float angleRandom, tiltRandom;
            float rand = Parameters.EmitterRate.Y * 0.15f;
            int x = (int)Math.Sqrt(Parameters.EmitterRate.Y);

            float degToRad = (float)(Math.PI / 180);

            for (int n = 0; n < Parameters.EmitterRate.Y; n++)
            {
                angleRandom = (float)((random.NextDouble() * 2.0 - 1.0) * Parameters.EmitterSpread.X);
                tiltRandom = (float)((random.NextDouble() * 2.0 - 1.0) * Parameters.EmitterSpread.Y);
                direction.X = (float)(Math.Cos((Parameters.EmitterAngle.X + angleRandom) * degToRad) *
                                      Math.Sin((Parameters.EmitterAngle.Y + tiltRandom) * degToRad) *
                                      Parameters.EmitterAngle.Z);
                direction.Y = 0;
                direction.Z = (float)(Math.Cos((Parameters.EmitterAngle.Y + tiltRandom) * degToRad) *
                                      Parameters.EmitterAngle.Z);

                position = Parameters.EmitterPosition;
                position.X += spacing * (n / x);
                position.Y = 0;
                particle = particles[AddPointReuse()];
                particle.Position = position;
                particle.Velocity = direction;
                particle.VelocityEvaluative = direction;
                particle.Age = 0;
                particle.Color = new Color3(time / 10.0f, time / 5.0f, time / 4.0f);
            }
        }

        public virtual void AddVolume(Vector3 min, Vector3 max, float spacing)
        {
            FluidParticle point;
            float dx = max.X - min.X;
            float dy = max.Y - min.Y;
            float dz = max.Z - min.Z;

            float randomDist = spacing * .1f;

            for (float z = min.Z; z <= max.Z; z += spacing)
            {
                for (float y = min.Y; y <= max.Y; y += spacing)
                {
                    for (float x = min.X; x <= max.X; x += spacing)
                    {
                        point = GetPoint(AddPointReuse());
                        point.Position = new Vector3(x + (float)random.NextDouble() * randomDist - randomDist / 2,
                                                     y,
                                                     z + (float)random.NextDouble() * randomDist - randomDist / 2);
                        point.Color = new Color3((x - min.X) / dx,
                                                 (y - min.Y) / dy,
                                                 (z - min.Z) / dz);
                    }
                }   
            }
        }


        // Ideal grid cell size (gs) = 2 * smoothing radius = 0.02*2 = 0.04
        // Ideal domain size = k*gs/d = k*0.02*2/0.005 = k*8 = {8, 16, 24, 32, 40, 48, ..}
        //    (k = number of cells, gs = cell size, d = simulation scale)
        public void Grid_Setup(Vector3 min, Vector3 max, float simScale, float cellSize, float border)
        {
            Vector3 borderVector = new Vector3(border, border, border);
            float world_cellsize = cellSize / simScale;
            
            GridMinimum = min;
            GridMinimum -= borderVector;
            GridMaximum = max;
            GridMaximum += borderVector;

            GridSize = GridMaximum;
            GridSize -= GridMinimum;
            gridCellSize = world_cellsize;
            
            // Determine grid resolution
            GridResolution = new Vector3((float)Math.Ceiling(GridSize.X / world_cellsize),
                                         (float)Math.Ceiling(GridSize.Y / world_cellsize),
                                         (float)Math.Ceiling(GridSize.Z / world_cellsize));
            
            // Adjust grid size to multiple of cell size
            GridSize = GridResolution * cellSize / simScale;

            GridDelta = new Vector3(GridResolution.X / GridSize.X,
                                    GridResolution.Y / GridSize.Y,
                                    GridResolution.Z / GridSize.Z);

            int numGridCells = (int)(GridSize.X * GridSize.Y * GridSize.Z);

            grid.Clear();
            gridCount.Clear();
            for (int i = 0; i < numGridCells; i++)
            {
                grid.Add(-1);
                gridCount.Add(0);
            }
        }

        public void Grid_Create()
        {

        }

        public void Grid_InsertParticles()
        {
            int gs;
            int gx, gy, gz;

            particles.ForEach(p => p.Next = -1);

            for (int i = 0; i < grid.Count; i++)
            {
                grid[i] = -1;
                gridCount[i] = 0;
            }

            int n = 0;
            foreach (var particle in particles)
            {
                //determine the grid cell
                gx = (int)((particle.Position.X - GridMinimum.X) * GridDelta.X);
                gy = (int)((particle.Position.Y - GridMinimum.Y) * GridDelta.Y);
                gz = (int)((particle.Position.Z - GridMinimum.Z) * GridDelta.Z);
                gs = (int)((gz * GridResolution.Y + gy)* GridResolution.X + gx);
                if (gs >= 0 && gs < grid.Count)
                {
                    particle.Next = grid[gs];
                    grid[gs] = n;
                    gridCount[gs]++;
                }
                n++;
            }
        }

        public void Grid_Draw(Matrix viewMatrix)
        {

        }

        public List<int> Grid_FindCells(Vector3 p, float radius)
        {
            List<int> gridCells = new List<int>();

            //if (p.X > Parameters.VolumeMax.X ||
            //    p.Y > Parameters.VolumeMax.Y ||
            //    p.Z > Parameters.VolumeMax.Z ||
            //    p.X < Parameters.VolumeMin.X ||
            //    p.Y < Parameters.VolumeMin.Y ||
            //    p.Z < Parameters.VolumeMin.Z)
            //{
            //    throw new InvalidOperationException("Point is out of bounds");
            //}

            for (int i = 0; i < 8; i++)
                gridCells.Add(-1);
            Vector3 sphereMin = new Vector3();

            sphereMin.X = (int)((-radius + p.X - GridMinimum.X) * GridDelta.X);
            sphereMin.Y = (int)((-radius + p.Y - GridMinimum.Y) * GridDelta.Y);
            sphereMin.Z = (int)((-radius + p.Z - GridMinimum.Z) * GridDelta.Z);

            if (sphereMin.X < 0)
                sphereMin.X = 0;
            if (sphereMin.Y < 0)
                sphereMin.Y = 0;
            if (sphereMin.Z < 0)
                sphereMin.Z = 0;

            gridCells[0] = (int)((sphereMin.Z * GridResolution.Y + sphereMin.Y) * GridResolution.X + sphereMin.X);
            gridCells[1] = gridCells[0] + 1;
            gridCells[2] = (int)(gridCells[0] + GridResolution.X);
            gridCells[3] = gridCells[2] + 1;

            if (sphereMin.Z + 1 < GridResolution.Z)
            {
                gridCells[4] = (int)(gridCells[0] + GridResolution.Y * GridResolution.X);
                gridCells[5] = gridCells[4] + 1;
                gridCells[6] = (int)(gridCells[4] + GridResolution.X);
                gridCells[7] = gridCells[6] + 1;
            }
            else
            {
                gridCells[4] = -1;
                gridCells[5] = -1;
                gridCells[6] = -1;
                gridCells[7] = -1;
            }
            if (sphereMin.X + 1 >= GridResolution.X)
            {
                gridCells[1] = -1;
                gridCells[3] = -1;
                gridCells[5] = -1;
                gridCells[7] = -1;
            }
            if (sphereMin.Y + 1 >= GridResolution.Y)
            {
                gridCells[2] = -1;
                gridCells[3] = -1;
                gridCells[6] = -1;
                gridCells[7] = -1;
            }

            return gridCells;
        }

        public int Grid_FindCell(Vector3 p)
        {
            int gc;
            Vector3 cell = new Vector3();
            cell.X = (int)(p.X - GridMinimum.X) * GridDelta.X;
            cell.Y = (int)(p.Y - GridMinimum.Y) * GridDelta.Y;
            cell.Z = (int)(p.Z - GridMinimum.Z) * GridDelta.Z;
            gc = (int)((cell.Z*GridResolution.Y + cell.Y)*GridResolution.X + cell.Z);
            if (gc < 0 | gc >= grid.Count) return -1;
            return gc;
        }

        public FluidParticle GetFirstGridParticle(int gc, ref int particleIndex)
        {
            currentParticle = grid[gc];
            if (currentParticle == -1)
                return null;
            particleIndex = currentParticle;
            return particles[currentParticle];
        }

        public FluidParticle GetNextGridParticle(ref int particleIndex)
        {
            FluidParticle ret = null;
            if (currentParticle != -1)
            {
                ret = particles[currentParticle];
                particleIndex = currentParticle;
                currentParticle++;
            }
            return ret;
        }

        public int GetGridCell(int x, int y, int z)
        {
            return (int)((z * GridResolution.Y + y) * GridResolution.X + x);
        }

        public ushort[] GetNeighborTable(int n)
        {
            return neighbor[n];
        }


        #endregion

        #region Private Methods

        #endregion
    }
}
